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ABSTRACT 

We use a set of Stackel potentials to obtain a local approximation for an effective 
third integral in axisymmetric systems. We present a study on the feasibility and 
effectiveness of this approach. We have applied it to three trial potentials of various 
flattenings, corresponding to nearly ellipsoidal, disky and boxy density isophotes. In 
all three good fit to the potential requires only a small set of Stackel potentials, 

and the associated Stackel third integral provides a very satisfactory, yet analytically 
simple, approximation to the trial potentials effective third integral. 

Key words: galaxies: structure, galaxies: elliptical and lenticular, galaxies: kinemat- 
ics and dynamics, methods: analytical, methods: numerical 



1 INTRODUCTION 

It is well-known that axisymmetric potentials do not in gen- 
eral allow a global and isolating third integral in addition 
to the energy and the component of the angular momentum 
along the symmetry axis. However, for many axisymmetric 
potentials, numerical experience shows that the majority of 
orbits are constrained by an effective third integral over a 
large number of orbital periods. As appears from modeling 
observed galaxies, a third integral is often required to de- 
scribe their dynamical structure {e.g. Binney et al.l990; De- 
jonghe et al.l996). Hence, it seems unescapable that distri- 
bution functions depending a priori on three integrals should 
be available for general models of elliptical galaxies. 

In a stellar dynamical context, the problem of approx- 
imating the effective third integral is logically connected to 
the philosophy behind the construction method of dynami- 
cal models for galaxies: obviously, numerical models can do 
with numerical, or even implicit third integrals, while ana- 
lytical models need analytical approximations. 

Amongst the numerical modelling methods, those based 
on Schwarzschild's approach (Schwarzschild 1979) are un- 
doubtly the most direct ones. In a given potential, a library 
of orbits is computed, while the time-averaged properties of 
the orbits are stored. The orbit library is supposed to sample 
integral space. A reproduction of the observables is sought 
by combination of orbits, populated with a non-negative 
number of stars {e.g. Rix et al. 1997; van der Marel et al. 
1998). This powerful and very general technique allows to 
construct general three-integral models, without any explicit 
reference to the third integral, since in principle orbits can 
be labelled using any phase-space point on them. Spectral 



analysis offers the possibility to simplify the expression for 
the orbits and to reduce storage requirements for the orbits 
(Binney & Spergel 1982; Papaphilippou & Laskar, 1996 & 
1998; Valluri & Merritt 1998; Carpintero & Aguilar 1998). 

However, for this type of modeling methods, the com- 
putational cost is relatively high. Moreover, smoothing is 
necessary because the singular boundaries of the orbital den- 
sities produce a fairly awkward sum of numerical functions, 
and the distribution function is a numerical function out 
of which the physical content may not be very easily ex- 
tracted. The implicit reference to a third integral, which 
makes Schwarzschild methods so flexible, is at the same time 
a handicap for extracting information about the role played 
by this integral. 

Analytical techniques would be more explicit. But no 
expression for a global third integral is known in general ax- 
isymmetric potentials - except for Stackel potentials (c/. de 
Zeeuw 1985; Dejonghe & de Zeeuw 1988). Therefore analyti- 
cal approximations have been built using perturbation meth- 
ods. Systems deviating moderately from spherical symmetry 
have been considered, and an approximate integral derived, 
which reduces to the total angular momentum in the spheri- 
cal limit (see Petrou 1983). Gerhard & Saha (1991) report on 
the use of three different perturbation techniques applied to 
a model that is initially spherical. (1) Methods based on the 
KAM theory, which have a validity that is often restricted to 
very small perturbations. In practice they find indeed that, 
whatever their order, these methods fail to track changes 
in the phase-space topology that occur when spherical sym- 
metry is broken by a significant amount. (2) The averaging 
method (see e.g. Verhulst 1979; de Zeeuw & Merritt 1983) 
in their example also fails to track the box-orbits (L^ = 0) 



© 0000 RAS 



2 De Bruyne, Leeuwin, & Dejonghe 



that emerge in a flattened system, and it gives at first-order 
only a rough approximation to the third integral. (3) Finally, 
a resonant method using Lie transforms does yield a good 
approximation to the third integral. It tracks the new orbital 
families emerging around resonances, and may in principle 
be carried on to high order. The analytical expressions are 
however rather intricate, specially for orders higher than 1. 

The last technique has been used at first-order by 
Dehnen & Gerhard (1993) to construct three integral oblate 
models for the perturbed isochrone sphere. An application 
to the boxy E3-E4 galaxy NGC 1600 has recently been pre- 
sented by Matthias & Gerhard (1999). 

However, one may need to model a galaxy whith a po- 
tential that is far from any known integrable potential, so 
that global perturbation techniques become ineffective. An 
alternative approach is then to derive a so-called 'partial 
integral'. The idea is that a third integral that can bo easily 
determined for certain groups of orbits is probably also a 
good approximation for similar orbits, de Zeeuw, Evans & 
Schwarzschild (1996) and Evans, Hafner & de Zeeuw (1997) 
derived such a partial integral for thin and near-thin tube 
orbits in scalo-frcc potentials. 

Yet another well known approach is the use of separable 
potentials. One is led to this idea because motion around 
the origin can be expanded in a Taylor scries and, since 
the lowest-order non-zero terms arc quadratic, can be de- 
scribed as perturbed harmonic motion. Moreover, van de 
Hulst (1962) has shown that motion around the origin can 
be seen as (unperturbed) motion in a Stackel potential. Ba- 
sically, one then uses a much enlarged set of reference in- 
tegrable potentials instead of merely quadratic ones. The 
Stackel potential is found by fitting term to term the Tay- 
lor series for the original potential (up to quartic terms) to 
the expansion for the Stackel potential. Local fits of simi- 
lar nature have been performed by de Zeeuw & Lynden-Bell 
(1985) and Kent & de Zeeuw (1991), who applied it to the 
solar neighbourhood, by expanding the effective potential in 
the vicinity of circular orbits. 

However, such a local fitting does not work for all or- 
bits. For instance, it may fail for orbits with largo radial 
extent. In many cases, though, an orbit can still be approx- 
imated by motion in a Stackel potential provided that the 
potential is chosen to be a good average approximation to 
the true potential in the region covered by that particular 
orbit (see Kent & de Zeeuw 1991). A report on the calcu- 
lation of such a good average and global approximation to 
the Galactic potential, together with some indications for 
the numerical implementation, can bo found in Dejonghe & 
de Zeeuw (1988) and Batsleer & Dejonghe (1994). 

In this paper, wo want to improve the generality of the 
application of this family of potentials. We propose to obtain 
a set of local approximations to the potential using Stackel 
potentials, each of which is an average approximation to the 
original potential in some region. Rather than expanding 
around an equilibrium point, which would make the extent 
of the region where the approximation holds difiicult to han- 
dle, we partition integral space and fit the potential in the 
corresponding regions of configuration space. We use the QP 
method (Dejonghe 1989) to adjust each local Stackel poten- 
tial. Section 2 describes this. Once the set of local Stackel 
potentials is obtained, we proceed by checking the quality of 
the orbits representation, and of the approximation for an 



effective third integral, where it exists. Obviously, we do not 
expect to describe the eventual stochastic motion this way. 
Moreover, we may not be able to reproduce all minor orbital 
families. Those questions are also addressed in section 2. In 
section 3 the actual fitting procedure is explained. Its per- 
formance is illustrated with three trial potentials, which cor- 
respond to meiss densities with resp. nearly ellipsoidal, boxy 
and disky isophotes. These trial potentials are presented in 
section 4. The results can be found in section 5. Section 6 
gives a discussion on the application of the method and the 
conclusions are given in section 7. 



2 THE DESIGN OF THE SET OF STACKEL 
POTENTIALS 

2.1 The principle 

Orbits with a given set of integrals {E, J), with E the (pos- 
itive) binding energy and J the component of the angular 
momentum along the rotation axis, fill a certain volume S in 
space (a meridional section of it is shown in figure lb). The 
intersection of this volume with a meridional plane is called 
the zero velocity curve (ZVC) . One proves easily that orbits 
with a set of integrals {E' , J) , with E' > E fill a volume that 
is completely embedded in the volume filled by orbits with 
the set of integrals {E,J). Similarly, orbits with E and J' 
larger than J, have a ZVC that lies inside the ZVC defined 
by (E, J). This means that all orbits in the shaded rectan- 
gle labeled TZ in integral space (see figure la) will remain 
interior to the shaded region indicated as <S in Fig lb. Any 
distribution function defined over TZ in integral space can 
therefore be associated with a potential that needs only be 
defined in S. Houcc, a subdivision of the {E, J)-plane into a 
number of rectangles TZ, is equivalent to subdividing space 
into a number of bounded domains S. For each of these do- 
mains, a Stackel potential can be determined that is locally 
a good approximation to the original potential. If in the 
end, the whole {E, J)-pIane is covered by a set of rectangles 
T^(E,j), also space will be covered by a set of domains S(^e,j), 
and the original potential will be completely fit by a set of 
locally fitted Stackel potentials. 

2.2 Practically 

In practice, the specification of a set of rectangles TZ in in- 
tegral space is equivalent to the specification of a grid, each 
gridpoint of which being the 'upper left corner' of TZ. A grid 
in E can be constructed by dividing the system into equal 
mass shells and computing the circular orbit energy for each 
shell radius. For every value of E, we consider an equidis- 
tant grid in J, with 7 values in the interval ]0, Jmax(i5)[, the 
upper bound corresponding to the circular orbit with that 
energy (see also figure 2). 

In Stackel potentials, orbits are determined by three in- 
tegrals of motion {E, Jjis). Following van der Marel et al. 
(1998), the different third integrals are parametrized with 
angles LUi, that correspond to the angle between the hori- 
zontal axis and the radius to the most external of the two 
points where the orbit intersects the ZVC (see figure 2). 
The point of the ZVC where I3 = /smax is found. This 
is the point where the thin tube orbit touches the ZVC. 
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Figure 1. Panel (a): A rectangle TZ. in integral space, of which 
the 'upper left corner ' corresponds to a gridpoint {E,J). Panel 
(b): All orbits with E' ^ E and J' ^ J will remain interior to the 
region S. In practice, the local Stiickcl potential is determined in 
the rectangle S' . 




Figure 2. For one value of E in the grid : zero velocity curves 
corresponding to 7 different values for J and points used to de- 
termine 7 values for I-^. 



The coordinates {R, z) of that point determine the angle 
Wmax = TT — Arctan(2:/(ii — Rc(E)). The uji are chosen lin- 
early between and ojmax, from which the values for I^iuji) 
are derived. 

As for the fit itself, it suffices to perform the Stackel fit 
to the potential associated with a rectangle TZ(e,j) only in 
the region interior to the velocity curve S(e..j)- In practice, 
we fit within the smallest rectangle S'(^^ jj which encom- 
passes S(^E,J)- 

The construction of a set of potentials that covers the 
complete system potential is a multi-step process. It is most 
efficient to start with a global fit to the potential, i.e. fitting 
a Stackel potential in the most extended domain S' , deter- 
mined, e.g., by the data. This domain defines the smallest 
value for E. The following steps improve the fit to the poten- 
tial by performing a number of local fits in smaller domains. 

The choice of domains where a fit will be performed 
is in principle arbitrairy, and depends only on how well 



the Stackel potential has to approximate the system poten- 
tial (within the limitations of what a regular potential can 
achieve) . 



2.3 When is a fit considered to be a good fit? 

The fit is checked by comparing orbits in both potentials. 
Orbits start from the ZVC defined by {E,J), at the point 
determined by the value of I3. The integration of orbits 
uses a fourth-order Runge-Kutta scheme with variable time- 
steps, proportional to the smallest of the radial and az- 
imuthal periods, estimated respectively as Tr ~ R/Vr and 
Tn ~ 2tvR'^/J. This ensures conservation of energy over 
100 Tn with a precision better than 10~^. 

One can consider various criteria for the comparison: 

• Surfaces of section: The Stackel potential should gen- 
erate orbits that are very similar to those in the original 
potential. A first inspection can be done by comparing sur- 
faces of section (SoS's) in both potentials. What we mainly 
expect to find this way is: 

- the proportion of non-regular orbits. If there is no 
additional integral of motion besides E and J, the orbits 
should fill uniformly the area bounded by the zero ve- 
locity curve on their SoS (Richstone 1982). If there is a 
third effective integral (effective meaning over the time 
the orbits are integrated) , orbits are regular, and appear 
as curves on the SoS. Experience has taught that in ax- 
isymmetric models, stochastic orbits exist only for very 
small values of J. In models with a core, as considered 
here, these orbits are due to the large excursions in R 
and z for small J, which allow resonances between the 
two degrees of freedom (Merritt 1999). 

- the minor families of orbits which may be present 
around some resonances in the original potential, that 
our Stackel potential does not generate. We call them 
unrecovered minor families. The thin tube orbit at fixed 
{E, J) determines the main orbital family. This family 
encircles the origin on the {z, Vz) SoS, or a point {Rth, 0) 
on the (R,vr) SoS, where Rth{E,J) is the point where 
the thin-tube orbit of given {E, J) intersects the equa- 
torial plane. Minor orbital families may be generated by 
other resonances than the ones for small J, and will show 
up on a SoS as 'islands' surrounding other points. 

• Orbital densities: 

Since the Stackel potential is aimed to be used for mod- 
eling, which is essentially assembling orbits to fit a given 
density, a good reproduction of the orbital densities is im- 
portant. 

These orbital densities ^{R, z; E, J, J3) are functions of 
{R,z) for each given orbit, defined for any set (£, J, Is). 
The orbital densities for both potentials can be compared 
by measuring the fraction of the mass that is located in the 
same place in both potentials. If the total mass of every or- 
bit is normalised to 1, and the cells surface is constant over 
the grid, the mass fraction correctly located (MC) can be 
calculated as 



MC = 1 - 5M = 1 - 



\{Vor[Rk,Zl) - Us[Rk,Zl))\/2. (1) 



In this, 5M is the mass fraction that is not located in the 



© 0000 RAS, MNRAS 000, 000-000 



4 De Bruyne, Leeuwin, & Dejonghe 



same cell in both potentials, Voi stand for the original density 
and vs stands for the density associated with the Stackel 
potential. Double counting is avoided by dividing the sum 
by 2. 

To calculate this, each orbit determined by a given set 
(-E, J, /a), is integrated over a large number of periods, 
and its orbital density i/{R, z; E, J, I3) (normalised to 1) 
is computed. This is done on a rectangular grid cover- 
ing [Rmin [E, J),Rmax{E, J)] X [0, Zmax{E, J)], by evaluating 
the time fraction spent in each cell. 

In principle the orbit should be integrated until the orbital 
density reaches quasi-stationarity. This we may define by de- 
manding that during a given time interval the value Vki in 
each ceU {Rk, zi) has varied by less than some small fraction. 
In practice, it takes an extremely long time to reach station- 
arity for orbits close to a resonance m : n with m, n large, or 
having a moderate degree of stochasticity - these two cases 
being difficult to distinguish (c/Binney & Tremaine 1987, 
p.176). 

Also modeling based on the Schwarzschild method has 
to deal with stationarity (see e.g. Schwarzschild 1993; Woz- 
niak & Pfenniger 1997). Since modeling aims at constructing 
equilibrium models, only orbits for which the orbital density 
reaches stationarity in limited time [i.e. less than a Hubble 
time) need in principle to be taken into account. Accord- 
ingly, the non-stationary orbital densities would not need to 
be precisely reproduced in the Stackel potential. 

• Conservation of I3: If locally a Stackel potential is a 
good approximation, the Stackel third integral I3 should be 
approximately constant along the regular orbits. Since the 
goal is to provide a good approximation for an effective third 
integral for the largest possible number of orbits, the varia- 
tion on I3 should be made as small as possible. If possible, 
fits on smaller domains will be performed if this improves the 
approximation for the integral. The approximation is eval- 
uated through the majcimal variation of I3 along the orbit 
with respect to its initial value, SI3/I3. 

On the other hand, one can argue that what really mat- 
ters is a nominal variation of I3, &I3/ HimoiiE, J) : for a 
given [E, J), how precise can our identification of orbits be 
when we estimate the third integral of motion by I3 within 
[0, /3max(^', J)X! A related question is how precise this iden- 
tification should be or, alternatively: how large is the toler- 
ated variation on I3I 

The Stackel potentials will be used as basis for a dynam- 
ical model. One of the aims of these models is to reveal the 
role of the third integral. Hence it is important to know 
that all structure caused by the dependence on this third 
integral can be traced. Consequently, the Stackel potential 
should be precise enough, so that the uncertainty on the re- 
sulting approximation for I3 is not larger than the scale of 
the smallest feature in I3 of the distribution function. For 
example, in the limit where the dynamical model would turn 
out to be essentially 2-integral, there would be no need to 
set requirements on the constancy of I3. 

This makes determining the maximum tolerance for errors 
in the fit to the potential an iterative process, which includes 
the modeling of the galaxy. First, a fit to the potential has 
to be obtained, and a model has to be constructed with 
that potential. The results of the model will then indicate 
whether the error on I3 was small enough or not. 

• Topology of the orbits: 



As pointed out by Gerhard & Saha (1991), approximate 
conservation of I3 along orbits alone does not ensure that I3 
provides a correct labelling for orbits (see also de Zeeuw & 
Merritt 1983). One still needs to check that the topology of 
orbital tori is similar to that of constant (i?, J, /s) surfaces: 
for then each orbit may be uniquely associated to one value 
of I3. This is done by comparing surfaces of section for orbits 
in the original potential, and I3 = est curves for given {E, J). 



3 HOW TO OBTAIN A STACKEL POTENTIAL 

The fitting procedure is derived from a deprojection method 
for triaxial systems presented in a paper by Mathieu & De- 
jonghe (1996). In that paper, a family of potential-density 
pairs is presented that can be used as building blocks for 
triaxial mass models. The spatial mass density and the po- 
tential can be both expressed in terms of the same basis 
functions. 

Using this method, it would technically be possible to 
obtain a Stackel approximation in two ways: (1) fitting the 
spatial mass density (2) fitting the potential. 

The first way is not ideally suited for performing local 
fits since the calculation of a local value for the potential 
requires an integration of the density over the entire volume 
of the galaxy. In the case of a fit in a restricted part of space, 
there is obviously no control on the behaviour of the mass 
density at larger distances, and this can have serious impli- 
cations for the potential obtained through integration. The 
correspondence between the forces generated by the original 
and the fitted potential, is likely to be a more useful indi- 
cation for the quality of the approximation than the corre- 
spondence between the spatial densities. For a global fit on 
the other hand, it does make sense to construct a Stackel 
potential through a fit on the spatial density. 

The fitting procedure is based on a Quadratic Program- 
ming method. The aim is to find a linear combination of 
basis functions that provide a good approximation to the 
original function. When fitting the potential, the basis func- 
tions are 



Fir) 



GM 



(2) 



(d + rpy 

The basis functions and their coefficients are chosen by min- 
imizing the variable 



^ ] Wkt [/original ('"fc, 2:^) ^ /calculated (^^fc , ^< )] ^ 



(3) 



with k and I an index covering the grid points in the domain 
5'. The weights Wki can be used to give different relative 
weights to the points in the grid, / stands for the function 
to fit. 

In an ellipsoidal coordinate system —7 ^ —a ^ A, 
with a and 7 negative, and the focal distance A = -y/la — 7] 
(de Zeeuw 1985), a Stackel potential can be written as 

V{\u)^g>.F{\)+g,F{u), (4) 

with F an arbitrary function (here corresponding to (|^) and 

A + Q 

9^ = (5) 
We now present the relevant formulas from Mathieu & 
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Dejonghe (1996), adapted for an axisymmetric system, that 
lead to an expression for the density in terms of the basic 
functions. 

The expression for the axisymmetric density is 

piKi^) ^ gli''W + gl^'it^) + '2gxg,.TijW,\], (6) 

where xpl^, A] is the first order divided difference of -0: 



X-u 



(7) 



and tp'{X) is the derivative tp[X, A]. 

The connection between •ipi''') F{t) (Hke in (^)) is 
given by 



2nGiP(r) = 2(r+7)F'(r)-F(r)+2l^ [F(r) - F{-a)] .(8) 

T + a 

The density can be expressed in terms of the basic func- 
tions by means of a third order divided difference: 



p{X,u) = H[X,u,X,u], 
with H{t) defined as: 
H(r) = {r + af^{T), 



(9) 



(10) 



with iP{t) given in IJ^. 

The divided difference of order n — 1 of a function G(r) 
is a function of the divided difference of order n — 2, and is 
given by 

1 G[ti,T3, . . . ,T„] - G[t2,T3, . . . ,Tn] 

G[ri, . . . ,r„J = . (11) 

n — T2 



4 THREE TEST POTENTIALS 

We build a smooth flattened model by combining a few 
spherical harmonics yj°(&) (Binney & Tremaine 1987), as 
follows: 



pir,e) = C[poir)Y°{e,(b) + f3p2{r)Y2°{e,(b) 



+ 



Po{r) 



p2{r) 



(1 + r2)2 

^2 



(1 +r2)3 
pi{r) = p2{r). 



(12) 

(13) 

(14) 
(15) 



/3 is a parameter which determines the flattening (/3 < 
for oblate systems). We set C = 2/{-k^/tt) for a total mass 
equal to 1. This simple choice for the model ensures p oc r^* 
and a fairly constant flattening. The resulting density for 
P = —0.3 is shown in flgure ^. It is nearly ellipsoidal, with 
axis ratio b/a ~ 0.8. We will in the following refer to the 
corresponding potential as ^'di. 
The corresponding potential is 

^eu{r,e) = -47rGC[V'o(r)yo°(e,0) + /3V'2(r)y2°(e,0) 



+ 



(16) 



where the 0; are related to the pi terms by equation (2-122) 
in Binney & Tremaine (1987). 

As explained in paragraph § 2, we use a grid in integral 
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Figure 3. The spatial density, with its axis ratio, and potential 
corresponding to 'ton, given by the sum of harmonics given in 
equation ( p^ ) and equation (p^), with /3 = —.3. 



space that defines the domains for the fit. The grid contains 
16 points in E (noted Ei), 7 points in J (noted Ji) and 7 
angles LJi. 

For each of the domains, a fit can be produced, based 
on the of equation if the focal distance A is set. The 
choice of A of the spheroidal coordinate system affects the 
quality of the fit. An average A is obtained by computing 
a number of orbits in the original potential, and evaluating 
the focal distance for which sections of the A = est ellipses 
most closely follow the inner and outer envelopes -Rmin(2;), 
Rmaxiz) of the orbits. Of course, there is more variation in 
A for the more extended domains 5'. 

We started with a global fit for the potential 'I'cii, i.e. 
a fit in the largest domain 5' (extending to about 180 in R 
and in z), corresponding to the orbits that have (i?i6,0). We 
take J — instead of Ji in order to avoid that orbits with 
smaller E should not be covered by the potential because 
their J < {Ji)Eie- 

This combination {Ei(i,0) is obviously the most unfa- 
vorable for the fit, and the variations in I3 are expected to 
set a standard for the worst case: the maximum (5/3)niax for 
all 7 orbits with different LJi is taken as the maximum toler- 
able variation. We demand that the variation of I3 along the 
orbits with lesser spatial extent (i.e. E' > Eie, J' > and 
the 7 uii, with (£', J') on the grid in integral space) is not 
worse than this. In this case, it turns out that the variation 
on /a is larger for orbits with (£15, Ji). Thus the next fit is 
done in the domain S' determined by (£15 ,0). 

After this fit, (5/3)max is the error on I3 for the orbits 
with (£15, Ji) and this new value is now used to decide on 
the next fit. In this way, fits are done in the domains cor- 
responding to (Eiafi), {Ei5,0), (£14,0), (£11,0), {E2,0) and 
(£1,0). 

If the correspondence between the orbits in a fitted and 
original potential happens to be very bad, it could be due to 
a poor correspondence between the force components of both 
potentials. Therefore, prior to the calculation of the orbits, 
it is useful to check the behaviour of the force components in 
the fitted and the original potential. Figure ^ shows the con- 
tours for the force components generated by ^'di (full lines) 
with the Stackel potential for the domain S'{Eii, 0) (dashed 
lines) in overlay. Since the force components in both poten- 
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Figure 4. The contours for the two force components gener- 
ated by ^£,11 (full lines) and the Stackel potential for the domain 
5'(ll,0){dashed lines). 
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Figure 5. The (i?, J)-grid for ^olh with different symbols for 
different potentials. We used a set of 6 Stackel potentialsto ap- 
proximate the original potential. 



tials do not appear too different, it is sensible to proceed 
with the comparison of integrated orbits. 

In figure q the points in the {E, J)-grid for ^l/eii are dis- 
played as an example of how a set of potentials can be built; 
different symbols indicate the different potentials used. As 
can be seen from the different symbols, some potentials of 
the set are used for a small number of orbits, while others 
are used for a large number of orbits. The fact that the po- 
tential fitted in the largest domain is only used for orbits 
with only one value for E (indicated by the pluses) suggests 
that this method can really improve the approximation. Also 
the orbits that remain close to the centre seem to require a 
potential fitted in a limited part of space. 

As a second trial case, a boxy density, with potential 
^box, is obtained from a combination of harmonics as de- 
scribed in ( p^ ) and (jl6|) with (3 = —0.5, and has an axis 
ratio h/ac^ 0.6. The contours are shown in figure 0. We also 





Figure 6. The spatial density, its axial ration and the potential 
corresponding to "fBoxi given by the sum of harmonics given in 
equation ( p^ ) and equation (^), with (3 = —.5. 





15 20 



Figure 7. The spatial density, its axial ratio and the potential 
for a Miyamoto-Nagai potential with b/a ^ 0.7. The contours for 
the spatial density are strongly disky. 



used a 16b X 7./ X T/g grid and approximated "I/box with a 
set of 8 Stackel potentials. 

Also a Miyamoto-Nagai (MN) potential with interme- 
diate flattening, b/a ~ 0.7, is taken as trial potential. MN 
models exhibit very strorig diskiness in the density contours, 
as can be seen in figure [TI This means that we are consid- 
ering a difficult case, that is actually on the verge of being 
unrealistic for elliptical galaxies. 

The grid has the same dimensions as the one described 
in the previous section. The MN potential is approximated 
with 8 potentials, the set of Stackel potentialsis constructed 
following the same strategy. 



5 PRESENTATION OF THE RESULTS 
5.1 Surfaces of section 

Typical surfaces of section display {R, vr) at 2 = 0, for the 
orbits having > when they cross the equatorial plane 
or display {z,v^) for vii = Q, for those orbits having vii> Q 
just after crossing. 

The SoS's immediately tell us that for the three poten- 
tials, all orbits of our phase-space grid have an effective third 
integral. Stochasticity, if present, was mild and would have 
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Figure 8. Surfaces of section for typical orbits in the MN poten- 
tial (small dots) and the Stackel potential (large dots). 
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Figure 9. Cumulative distribution for the fraction of mass 'cor- 
rectly located' in the orbital densities. The full line is for 'I'cii, the 
dotted line is for ^box a^nd the dashed line for MN. 



required very long integrations to be visible on the SoS's. 
We did not carry out such integrations, because the SoS's 
were only meant as a preliminary check. In practice, any- 
way, mildly stochastic orbits behave very much like regular 
orbits, and we can hope to approximate them by an orbit in 
a Stackel potential. 

For tpaii and t/'box, all orbits in our library have similar 
SoS's in both potentials. There is no evidence for stochas- 
ticity, nor for unrecovered minor families. For the MN po- 
tential, all orbits also appear as regular, but we found a few 
unrecovered minor orbital families. 

Most of the orbits trapped around resonant tube or- 
bits - what we called 'minor orbital families' in §2.3 - were 
present in the local Stackel potential. 

In figure ^ SoS's of a few orbits are displayed for the MN 
potential (small dots) and the Stackel potential (large dots). 
Both panels show a 'minor resonance' evidenced by small 
'islands'. In the left panel, showing orbits with large I3, the 
orbits trapped around minor resonances are well reproduced 
within the Stackel potential. The right panel displays orbits 
with smaller I3. Two orbits are trapped around a 1:1 res- 
onance which does not exist in the Stackel potential. They 
are examples of unrecovered minor orbital families. They 
are always found for small values of I3. Therefore, they cor- 
respond to orbits remaining close to the equatorial plane, 
where the MN potential is very much distorted by a strong 
diskiness. 

Even for such extreme 'diskiness', these orbits only rep- 
resent < 1% of the orbit library for MN. These obviously 
are orbits that we will not be able to take into account when 
building a model that uses Stackel potentials. Although the 
unrecovered minor families are only a small fraction of the 
entire orbital library, we will exclude them from the statis- 
tics of the other checks discussed hereafter. As such, we are 
not taking those orbital families into account that no one 
expects to be approximated by a Stackel potential. 



5.2 Orbital weights 

For every orbit, the conserved mass fraction MC is calcu- 
lated following equation (|l]), using a 20 x 20 grid in {R,z). 
We integrate each orbit for at least 50 Tq , at most 200 Tq . 
We check every 50 Tn whether quasi-stationarity has been 
reached. Quasi-stationarity is assumed when for every cell 
the orbital density has varied by less than 1% in 50 Tn. 

The results are presented in figure |^ which gives the cu- 
mulative distribution for MC, the fraction of mass 'correctly 
located' in the orbital densities. A few orbital densities for 
MN have MC < 60%, these correspond to minor resonances 
that were not fitted by our Stackel potentials. Excluding 
those, the average MC is 99% for 98% for ^box and 

95% for MN. 



5.3 How constant is the Stackel /a? 

For every potential, the conservation of I3 along the or- 
bits is estimated by computing the relative variation 5^ = 
5l3{E,J,Is)/Is{E,J,Is). At given {E,J), the maximum of 
this quantity is derived among orbits with different I3: 

5n,ax = maX,= l,7[<5/3(-E', J, h,i)/h{E, J, /3,i)]. 

Figure |lo| shows how log((5J„ax) is affected as more po- 
tentials are added to the set approximating t/joII- 

The black dots represent the E ( J = 0) of the domain 
S' where the potentials are fitted. Dark gray shades repre- 
sent small variations on /a , light gray shades represent larger 
variations. For each new local potential that is considered, 
a number of previously void cells of the {E, J) grid are filled 
and a number of values are replaced by new ones. It is clear 
that adding new potentials to the set improves the conser- 
vation of the third integral along the orbits, the difference 
between the global fit (upper left panel of figure |lo|) and the 
results of the complete set (lower right panel) is remarkable. 

In the left panels of figure U/U we present histograms 
of 5"" for from top to bottom: ipaii, V'box and MN. For the 
MN potential, the resonances absent from the Stackel ap- 
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Figure 11. Histograms of 5'" in tlie left panels, <5" (i.e. 5/3 nor- 
malised to /3max(^'i •^)) in the middle panels and of SI3 nor- 
malised to ^SmaxC^) ™ the right panel. 



proximation are not considered, a complete orbital library 
contains 784 orbits. 

The average for is 2% for tpcu, while 75% of the orbits 
have < 1.7%. For i/^box the average for is ~ 2.6% and 
75% of the orbits have 5'' < 2.7%. For MN, the average 
value is ~ 10%, and 75% of the orbits have .5'' < 12%. The 



variation of 5^ along each orbit is of the order of the fitting 
error on the derivatives of the potential. 

In the middle panels of figure |l^ the results in terms 
of a nominal variation 5" = [8Li{E, J, 1^) / l2,nia.^{E , J)\ are 
shown. For the nearly ellipsoidal potential ^'di, 75% of the 
orbits have 5" < .3%, with an average value of ~ .2%. For 
^'box the average is ~ .7% while 75% of the orbits have 
5" < 1%. For MN, the average value is ~ 1.7% and 75% of 
the orbits have a nominal error on I3 less than 2.3%. 

Dehnen & Gerhard (1993) used a similar indicator, but 
they normalised SI3 to /3niax(£') , which gives smaller values 
for the nominal variation, as can be seen in the right panels 
of figure |l^. The average for respectively ^eii, ^'box and MN 
is .1%, .5% and 1.3%. 



5.4 Topology of constant I3 varieties 

We plot on the same graph (figure [l2| ): (i) the {R,R) SoS 
of orbits as they cross the z = plane (with i > 0) in the 
original potential; and (ii) the invariant curves defined by 
I3 — est (at z — 0). This is shown here for some intermediate 
value of E {E4 of our grid) and J = 0. Radial orbits are in 
principle difficult to map, because of the transition between 
loop and box orbits (see Gerhard & Saha 1991). They are 
however very well described here, both for "ifcu and ^l/box- 
The agreement is somewhat poorer in the case of MN for 
the transition region between box and loop orbits; however, 
the topology is still correctly described, in the sense that we 
can establish a one-to-one correspondence between the two 
sets of curves. This shows that /a may be used to label the 
orbits. 
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Figure 12. Comparison of SoS's and J3 level curves at given 
{E4, J = 0), for 'I'cii (upper panel), 'ifhox (middle panel), and 
MA'' (lower panel). Dots: surfaces of section (R, R) at z = with 
i > 0. Lines; the intersection of I3 = est surfaces with the z = 
plane. 

6 DISCUSSION 

For the construction of dynamical models, the use of Stackel 
potentials yields a number of advantages, amongst them (1) 
the property that the density and dynamical moments can 
be calculated analytically and (2) the absence of regularisa- 
tion problems. Finally, once the dynamical model has been 
completed, the distribution function is known in an easy to 
use and analytical form. Though these models require a spe- 
cial form of the potential, they offer a great flexibility during 
the further modeling process. In a Stackel potential there is 
a function of 1 variable that can be freely chosen. This free- 
dom is advantageous, and can, of course, be exploited at the 
fullest when performing a fit to the galaxy potential. 

Analogous to local Stackel potentials , dynamical mod- 
els built within this approximation will yield local distribu- 
tion functions, simply because we have adapted the coordi- 



nate systems in each domain S where a local Stackel poten- 
tialis constructed. The use of a QP-method (Dejonghe, 1989) 
for dynamical modeling gives a large freedom for the basis 
functions. In this case, we can use basis functions for the 
distribution function that contribute only in limited parts 
of integral space. In practice, these limited parts will corre- 
spond to the rectangles TZ from the local potentials. 

The fact that the space of Stackel potentials has mea- 
sure zero in the space of all axisymmetric potentials, raises 
the obvious concern whether Stackel potentials are suitable 
for the representation of any galaxy potential. 

The scientific objections cluster around the following 
arguments: 

(1) 'What about central cusps?' Indeed, central density 
cusps generally cannot be generated by a Stackel potential in 
an ellipsoidal coordinate system. However, Sridhar & Touma 
(1997) showed that a Stackel potential in parabolic coordi- 
nates can create a central cusp. Whether these potentials 
can be used in the scheme we propose here remains to be 
investigated. On the other hand, it may well be possible 
that there is no point in providing a good approximation 
to the third integral in regions where the cuspiness of the 
density becomes important. Indeed, it is likely that the cen- 
tral density cusps are related to the presence of a central 
massive object (see e.g. review by Kormendy & Richstonc 
1995; Richstone et al. 1998). This may cause a higher frac- 
tion of ergodic orbits over a Hubble time. Such orbits will fill 
regions of phase-space where the distribution function de- 
pends essentially on the energy, or, for less ergodic regions, 
also on angular momentum (c/Merritt 1999). The fraction 
of regular to ergodic motion in a real galaxy stands as a ma- 
jor unknown. The most popular scenario at the moment is 
that ergodicity has gradually been introduced in the system, 
as orbits that pass close to the central mass concentration 
were perturbed by it (Gerhard & Binncy 1985; Valluri & 
Merritt 1998). One may then suppose that the original or- 
bital tori which have not been disrupted still underlie most 
of the features in phase-space- with essentially homogeneous 
6-D density between them. 

(2) 'Docs a Stackel approximation allow an accurate fit 
to a typical galactic potential?' As was already noticed for 
some time, the main orbit families found by numerical inte- 
gration in general triaxial potentials are present in a Stackel 
potential (Schwarzschild 1979; de Zeeuw 1985), but obvi- 
ously there is no place in an integrable potential for smaller 
orbital families nor stochastic orbits. These minor orbital 
families appear to occupy only a small volume fraction of 
phase-space, as long as figure rotation is unimportant (Ger- 
hard, 1985; Binncy, 1987) In practice, Stackel potentials do 
turn out to provide reasonably good global fits for systems 
without central mass concentration (Dejonghe et al. 1996) 
or for regions beyond the influence of the central mass con- 
centration (Emscllcm et al. 1999). It is true that the use of 
one single Stackel potential assumes a single confocal sys- 
tem where the streamlines of the mean stellar motion coin- 
cide with the coordinate lines of the ellipsoidal coordinate 
system. Hence, the construction of such a single Stackel po- 
tential inevitably involves some sort of averaging (dc Zeeuw 
& Lynden-Bell 1985; Dejonghe & de Zeeuw 1988), and is 
not always considered to be a sufiiciently general approach 
(at least in triaxial cases, Binney 1987; Merritt & Fridman 
1996). However, the use of confocal streamlines seems to be 
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a valid assumption for individual orbits (Anderson & Statler 
1998). Moreover, studying logarithmic potentials, those au- 
thors find that the mean velocities can be well fitted using 
local coordinate systems with confocal coordinate lines, with 
the focal distances taken within a fairly narrow range. Our 
approach allows to carry this idea one step further, on the 
level of the distribution function. 



7 CONCLUSIONS 

In this paper, we present an extension of the available ap- 
proximations for the effective third integral in axisymmetric 
systems, obtained by using a set of Stackel potentials as 
representation for a galaxy potential, instead of one single 
Stackel potential. We have studied the feasibility and effec- 
tiveness of this method. 

The creation of a set of local Stackel potentials is done 
through a fit on the system potential. There is a large free- 
dom in the choice of the domains where the local fits are 
done and the composition of the set of Stackel potentials. 
The set of approximating local potentials can be extended 
until the desired precision is obtained. 

We have tested the method on three potentials: (A) 
a harmonical potential (^'oii) that behaves very smoothly 
with a nearly ellipsoidal density; (B) a harmonical poten- 
tial (^box) with a boxy density; and (C) a Miyamoto-Nagai 
model (MN) with a density that has an exaggerated disky 
structure. 

As one may expect, the model with smallest diski- 
ness/boxiness is easiest to approximate using sets of Stackel 
potentials. But what really motivates this study is to check 
the estimate made for the effective third integral, if we ap- 
proximate it by the third integral I3 that the Stackel poten- 
tials define. 

The quality of the approximation is checked in several 
ways: (1) surfaces of section, that reveal possible resonances 
and irregular orbits, (2) conservation of orbital weights, 
which is important for the reconstruction of the spatial den- 
sity, (3) conservation of /a along the orbits, in order to vali- 
date the labeling of orbits, (4) topology of the orbital space. 
The conservation of orbital weights and the conservation of 
I3 are the criteria that can be best expressed numerically. 

According to the expectations, the approximation of 
^'cii has orbits that are very similar to the orbits in the 
original potential, as can be judged from the surfaces of 
section and the topology. The results for the conserva- 
tion of orbital weights (an average of 99% for 'I'cii) and 
I3 { S"- = Sh{E,J,h)lh{E,J,h) ~ 2% and 5" = 
[Sh{E, J, /3)//3max(S, J)] ~ .2% for *eii) coufirm the qual- 
ity of the approximation. 

The method also proves to be successful for V&box, with 
an average of 98% for the conservation of orbital weights, 
6'' ~ 2.6% and S" ~ .7%. The SoS's and the topology of the 
orbits also confirms this good result. 

The success of the method on the Miyamoto-Nagai po- 
tential is somewhat less, because of the strong diskiness. 
Still, using a small set of Stackel potentials, we are able to 
reproduce most orbits with satisfactory accuracy, except for 
a few resonances. The resonances that are not reproduced 
by the set of Stackel potentials, all have small values of ^3, 
i.e. the orbits lie in the region where the diskiness is im- 



portant. Leaving this resonances out of consideration, the 
orbitals weights seem to be well conserved (on average 95% 
for MN) and the topology of the orbits is well reproduced. 
Also for these potentials, the Stackel I3 can be used as label 
for the orbits (5'" ~ 10% and 5" ~ 1.7%). The potentials of 
observed elliptical galaxies are generally rounder than the 
strongly disky Miyamoto-Nagai models. 

Given the positive results found for the three rather 
different models considered, it seems to be possible, using 
a reasonable number of local Stackel potentials, to provide 
good approximations for a third integral, suitable for label- 
ing orbits in dynamical models. 

With respect to existing approximations for a third in- 
tegral in axisymmetric systems, the advantage of this new 
one is mainly that, while it can be envisaged for application 
to systems with arbitrary flattening, it also yields a sim- 
ple analytic expression for the approximate third integral 
locally. This method will be fully exploited if it is used to 
build, for roughly axisymmetric galaxies, explicit distribu- 
tion functions that depend on three integrals. 



REFERENCES 



Anderson R.F., Statler T.S., 1998, ApJ, 496, 706 
Batsleer P., Dejonghe H., 1994, A&A, 287, 43 
Binney J. J., 1987, in Structure and Dynamics of Elliptical Galax- 
ies, T. de Zeeuw, ed., (Dordrecht: Reidel), p229 
Binney J. J., Spergel D., 1982, ApJ, 252, 308 

Binney J. J., Tremaine S., 1987, Galactic Dynamics, Princeton 

University Press, Princeton 
Binney J. J., Davies R.L., lUingworth G.D., 1990, ApJ, 361, 78 
Carpintero D.D., Aguilar L.A., 1998, MNRAS, 298, 1 
Dehnen W., Gerhard O.E., 1993, MNRAS, 261, 311 
Dejonghe H., 1989, ApJ, 343, 113 
Dejonghe H., de Zeeuw, P.T., 1988, ApJ, 329, 720 
Dejonghe H., De Bruyne V., Vauterin P., Zeilinger, W.W., 1996, 

A&A, 306, 363 
de Zeeuw P.T., 1985, MNRAS 216, 273 
de Zeeuw P.T., Lynden-Bell D., 1985, MNRAS, 215, 713 
Merritt D., 1983, ApJ, 267, 571 
, Evans N.W., Schwarzschild M., 1996, MNRAS, 



MNRAS, 303, 495 
1997, MNRAS, 286, 



de Zeeuw P.T. 
de Zeeuw P.T 
280, 903 

Emsellem E., Dejonghe H., Bacon R., 1999, 
Evans N.W., Hiifner R.M., de Zeeuw P.T., 
315 

Gerhard O., 1985, A&A, 151, 279 

Gerhard O. , Binney J., 1985, MNRAS, 216, 467 

Gerhard O.E., Saha P., 1991, MNRAS, 251, 449 

Kent S.M., de Zeeuw P.T., 1991, AJ, 102, 1994 

Kormendy, J., Richstone, D. 1995, ARA& A 33, 581-624. 

Mathieu A., Dejonghe H., 1996, A&A, 314, 25 

Matthias M., Gerhard O., MNRAS, submitted, istro-ph/9901036 

Merritt D., 1999, PASP, 111, 129 

Merritt D., Fridman T., 1996, ApJ, 460, 136 

Papaphilippou Y., Laskar J., 1996, A&A, 307, 427 

Papaphilippou Y., Laskar J., 1998, A&A, 329, 451 

Petrou M., 1983, MNRAS, 202, 1195 

Richstone D.O., 1982, ApJ, 252, 496 

Richstone D.O., Ajhar E.A., Bender R., Bower G., Dressier A., 
Faber S.M., Filippenko A.V., Gebhardt K., Green R., Ho L.C., 
Kormendy J., Lauer T.R., Magorrian J., Tremaine S., 1998, 
Nature 395, 14 

Rix H.-W., de Zeeuw P.T., Cretton N., van der Marel R.P., Car- 
ollo CM, 1997, ApJ, 488, 702 



© 0000 RAS, MNRAS 000, 000-000 



Approximate third integrals for axisymmetric potentials using local Stdckel fits 



Schwarszchild M., 1979, ApJ 232, 236 

Schwarszchild M., 1993, ApJ, 409, 563 

Sridhar S., Touma J., 1997, MNRAS 292, 657 

Valluri M., Merritt D., 1998, ApJ, 506, 686 

van de Hulst H.C., 1962, Bull. Astr. Inst. Neth., 16, 235 

Verhulst F, 1979, RSLPT, 290, 435 

van dcr Marcl R.P., Crctton N., dc Zccuw T.P., Rix H.-W., 1998, 

ApJ, 493, 613 
Wozniak H., Pfenniger, D., 1997, A&A, 317, 14 



© 0000 RAS, MNRAS 000, 000-000 



